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Modeling complex chemical effects in 
turbulent nonpremixed combustion 

By Nigel S. A. Smith 


1. Motivation 

Virtually all of the energy derived from the consumption of combustibles occurs in 
systems which utilize turbulent fluid motion. Since combustion is largely related to 
the mixing of fluids and mixing processes are orders of magnitude more rapid when 
enhanced by turbulent motion, efficiency criteria dictate that chemically powered 
devices necessarily involve fluid turbulence. 

Where combustion occurs concurrently with mixing at an interface between two 
reactive fluid bodies, this mode of combustion is called nonpremixed combustion. 
This is distinct from premixed combustion where flame-fronts propagate into a ho- 
mogeneous mixture of reactants. These two modes are limiting cases in the range of 
temporal lag between mixing of reactants and the onset of reaction. Nonpremixed 
combustion occurs where this lag tends to zero, while premixed combustion occurs 
where this lag tends to infinity. Many combustion processes are hybrids of these 
two extremes with finite non-zero lag times. 

Turbulent nonpremixed combustion is important from a practical standpoint be- 
cause it occurs in gas fired boilers, furnaces, waste incinerators, diesel engines, gas 
turbine combustors, and afterburners etc. To a large extent, past development of 
these practical systems involved an empirical methodology. Presently, efficiency 
standards and emission regulations are being further tightened (Correa 1993), and 
empiricism has had to give way to more fundamental research in order to understand 
and effectively model practical combustion processes (Pope 1991). 

A key element in effective modeling of turbulent combustion is making use of a 
sufficiently detailed chemical kinetic mechanism. The prediction of pollutant emis- 
sion such as oxides of nitrogen (NO x ) and sulphur (50 x )? unburned hydrocarbons, 
and particulates demands the use of detailed chemical mechanisms. It is essential 
that practical models for turbulent nonpremixed combustion are capable of handling 
large numbers of ‘stiff’ chemical species equations. 

1.1 Reactive Species Closure problem 

A common way of idealizing a turbulent flow field is to decompose it into an 
averaged flow component and a deviational contribution. The nature of the devia- 
tional component depends upon the flow and the averaging scheme, but the object 
of the decomposition is to be able to understand and predict the development of 
the average component without detailed knowledge of the deviations present in each 
realization of the flow. 

The classical difficulty faced in modeling turbulent nonpremixed combustion is 
that of closing the averaged equations for chemically reactive species. The instan- 
taneous equation for the evolution of the mass fraction Y a of a reactive species a is 
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the following, 


+ £r(' m ' Y ° ) = £- ipD °i£i ) + ** • (i) 

where w Q is the net chemical production rate of the species a, and D Q is the corre- 
sponding molecular diffusivity where a simplified Fickian approximation has been 
made to model molecular transport. Applying a traditional averaging scheme, such 
as density weighted (Favre) unconditioned ensemble averaging, yields the following, 

In order to close the averaged species equation a model must be provided for 
the averaged source term w Q . First order closures that evaluate the instantaneous 
chemical rate expressions with averaged species concentrations and temperature, 

w q (y i , y 2 , . . . , y n , t) « u> Q (yi ,y 2 ,...,y n ,t) (4) 

are known to be highly inaccurate in combustion cases of practical interest. The 
chemical reactions encountered in combustion processes are highly nonlinear, and 
thus small perturbations in the input parameters can cause very large changes in 
the computed reaction rate. 

1.2 Conditional Moment Closure method 

The philosophy underpinning the Conditional Moment Closure (CMC) method, 
as described by Bilger (1991, 1993), is to minimize the level of perturbations from 
the mean by averaging the reactive species equations conditionally upon a conserved 
scalar mass fraction. In so doing, the resultant statistical moments account for the 
variations in fluid concentration which result from turbulent mixing alone. At the 
expense of adding an additional computational dimension to the modeling prob- 
lem, conditional averaging allows chemical closure to be achieved for most cases of 
nonpremixed turbulent combustion. 

The average of a fluctuating turbulent quantity A, conditional upon the conserved 
scalar mixture fraction £(x;,tf) being equal to a sample value 77 , is the following (see 
Klimenko 1990): 


< A(xi,t) I Z(xi,t) = T) >= d- J J A(xi,t)6({(xi,t) - T))dxidt (5) 

In the above definition, P v is the probability density function of the conserved 
scalar at the location x, and time t, and 6 denotes the Dirac delta function. In all 
that follows, the full conditional averaging operator < ... | £(x,, t) = tj > will be 
abbreviated to < . . . | r? > for the sake of brevity. 
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Klimenko (1990, 1992) and Bilger (1991, 1993) independently showed that the 
evolution of the conditional mean mass fraction Q a =< Y Q \ r\ > of a reactive 
species a is governed by the following, 


, dQ a , ^ ^ dQ a 1 ^ ^ d 2 Q a 


f < pw Q | 1 > +e q (6) 


where the molecular diffusivities of all species are assumed to be uniform. The 
residual term e q is a conditional correlation between deviational velocity and mass 
fraction, which is typically assumed to be small. This assumption is not valid in 
cases where substantial premixing of the reactants occurs, in which case this term 
can become very important in the reactive species equation (Bilger 1991). 


e? = ~Q^( P v <P\v>< uWa I V >)/ P V 


( 7 ) 


The symbol x denotes the instantaneous scalar dissipation rate and is defined 
(below) in terms of the mixture fraction £. 




(S) 


In order to close the CMC scalar equation, means of determining < x I V > an d 
< Wa \ V > are required. 

Klimenko (1992) (see also Klimenko and Bilger 1992) showed that the conditional 
mean scalar dissipation rate should be determined from the conserved scalar PDF 
equation to ensure conservation of mass. 

d d 1 ^ 

gj(< P I V > P v) + ^"(< PUi\n> Pv) = PX ' 77 > P ^ + 

The residual term describes molecular diffusion of the PDF in physical space and 
is negligible at high Reynolds numbers. 

Closure of the chemical source term is achieved via a simple first order ap- 
proximation involving conditionally averaged species mass fractions < w a \ r) 
Wa(Qu • • • 5 Qn, < T | 77 >). This closure approximation is valid in all cases except 
where the combustion system is close to extinction, since in those cases deviations 
from the conditional means are large. In such instances doubly conditioned mod- 
eling, using conditions upon mixture fraction and a reaction progress variable, is a 
suitable course of action (Bilger 1991). 

The single greatest advantage of the CMC method over other turbulent non- 
premixed combustion models is its ability to cope with very detailed chemical de- 
scriptions. When contrasted with the Joint Probability Density Function (JPDF) 
method (see Pope 1985, 1991) the CMC method only adds a new equation to be 
solved with each new species rather than an additional dimension to the problem. 
In comparison to laminar flamelet methods, the CMC method does not require 
there to be a large separation of characteristic scales between mixing and chemical 
reaction in order to be applicable. 
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2. Objectives 

The objective of this study was to evaluate the CMC method for turbulent non- 
premixed combustion, verify or debunk current model assumptions, and where pos- 
sible suggest model refinements. 

Direct numerical simulations (DNS) of nonpremixed combustion in isotropic de- 
caying turbulence were carried out as a first step in this study. By studying non- 
premixed combustion in isotropic turbulence with an isotropic flame distribution, 
all of the computational grid points can be used in calculating spatially degener- 
ate statistics. These statistics allow important modeling problems to be examined 
without interference from the extraneous, though also important, complications 
introduced by mean gradients. 

Of key interest in this study is the effectiveness of the first order chemical closure 
employed by the current CMC methodology in a system which embodies the es- 
sential elements of combustion chemistry. In this study one- and two-step chemical 
mechanisms were employed to describe H 2 — N 2 fuel burning in air. The one-step 
mechanism is composed of the non-carbon step from the two-step wet-CO mech- 
anism of Chung and Williams (1990), and allows for global reaction termination 
due to radical consumption even though no radical species are actually carried (see 
below). 

2H 2 + 0 2 ->2 H 2 0 (/') 

The two-step mechanism carries the crucial radical, monatomic hydrogen ( H ), 
and contains distinct reactions for radical formation and consumption. The chemical 
rate constants for reactions I and II (see below) were also taken from the non-carbon 
steps of Chung and Williams (1990). 

3H 2 + 0 2 ^2H 2 0 + 2H (I) 

H + H + M->H 2 +M (II) 

In addition to the abundance of evidence in chemical kinetic literature, recent DNS 
studies conducted at the CTR have demonstrated the importance of the explicit 
calculation of chemical radical species in combustion simulations (see Mantel 1994, 
Vervisch 1992). 

Differential diffusion of reactive scalars in nonpremixed combustion is left unmod- 
eled by both the CMC (Smith 1994) and JPDF methods (Yeung and Pope 1993). 
Although it is of diminished importance in highly turbulent combustion, practical 
cases have been found where differential diffusion is significant (see Chen et al 1990, 
Smith et al 1993, Smith 1994, Bilger 1982). The secondary objective of this study 
was to observe differential diffusion phenomena by comparing simulations that are 
identical save for the specification of uniform or non-uniform species molecular dif- 
fusivities. These observations are to be used to develop model refinements for the 
prediction of differential diffusion behavior. 
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2.1 Modeling method 

The spatially degenerate CMC and PDF equations corresponding to statistically 
isotropic conditions with uniform molecular diffusivities are given below, 


dQ a 

dt 


2 < X I V > 


d 2 Q a 

a /? 2 


+ < W a I T) > 


( 10 ) 


^(< P \ V > Pv) = PX \ V > Pr,) ( 11 ) 

where the conditional averages are taken over the entire domain and residual terms 
have been neglected. 

The CMC equation (Eq. 10) was solved with the conditional mean scalar dissi- 
pation rate profile being given by the PDF equation in the following manner: 

< * 1 ” >= < „ |~, 2 > p, // M (< 1 5 (12) 


In this study it is possible to use PDF information from the simulation to determine 
the conditional mean scalar dissipation rate, but in practice this information will 
not typically be available to the modeler. Following the approximations used in 
practice, the conserved scalar PDFs were assumed to be clipped Gaussian in form. 
This assumption reduces the number of degrees of freedom in the PDF to two, 
namely specification of the conserved scalar mean and variance. Beta functions are 
an arguably superior assumed form (Girimaji 1991), however Gaussians readily lend 
themselves to accurate integration (Smith 1994). 

In cases where differential diffusion is significant, further terms appear in the 
CMC and PDF equations. Discussion of differential diffusion modeling is deferred 
to Section 4. 


2.2 Simulation conditions 

An upgraded version of the code used by Ruetsch (1994) was employed in the 
direct numerical simulations. The new code includes a flexible multi-step chemical 
kinetics module for handling arbitrary thermochemistry. The code retains original 
features such as the high-order compact finite differencing scheme described by 
Lele (1992) for spatial differencing, and the third order Runge-Kutta timestepping 
algorithm of Wray. The Navier-Stokes characteristic boundary conditions described 
by Poinsot and Lele (1992) are also retained, but this feature was not used in this 
study due to the periodic nature of the simulation domain. 

The simulations performed to date have been two dimensional (129x129) in order 
to extensively test the new code and simulation conditions, before expending a great 
deal of computation time on three dimensional simulations. 

The turbulent field was initialized using an incompressible phase scrambled ki- 
netic energy spectrum for the velocity components and a conserved scalar. The 
initialized conserved scalar field can be seen Fig. 1 where black regions denote pure 
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FIGURE 1. Initial distribution of the conserved scalar. White regions denote f = 1 
while black regions denote £ = 0. Initial conserved scalar unmixedness ft = 0.84. 


oxidizer zones and white regions denote pure fuel zones. Scalar unmixedness can be 
defined as, 

Q=<e> /(< t > (i- < t >)) (i4) 

which can be seen to be a normalized measure of the fluctuation level. Unmixedness 
varies between zero, where the scalar field is homogeneous, and unity where only 
pure fuel and pure oxidizer zones exist with no mixing at all. The initial conserved 
scalar fields used here had initial unmixednesses of ft « 0.8 in all cases. 

Reactive species mass fractions and internal energy were mapped onto the con- 
served scalar field using adiabatic chemical equilibrium relationships between mix- 
ture fraction (conserved scalar) and the reactive scalars. The adiabatic equilibrium 
reactive scalar mass fraction profiles are plotted versus conserved scalar mixture 
fraction in Fig. 2. Note that the fuel is comprised of almost 97% nitrogen (A^) by 
mass, thereby giving a stoichiometric mixture fraction of £ s toic — 0.5. As the mean 
mixture fraction for the simulations was also 0.5, the overall equivalence ratio was 
unity in all cases. 

The imposition of a 4 hot’ scalar field onto the initial ‘cold’ field solution required 
the adjustment of the density field to minimize the effect on the pressure field. Al- 
though the scalar mapping was essentially a constant pressure process, the resultant 
pressure field had a noticeable acoustic component that arose from the imposed im- 
balance in the momentum equations. This acoustic adjustment led to an initial 
root mean squaxe pressure fluctuation that was approximately 5.e — 5 of the mean 
pressure. Due to periodicity, the acoustic waves were unable to leave the domain 
but did slowly decrease with time due to dissipation. 
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FIGURE 2. Initial conditional mean reactive scalar profiles for two-step chemical 
description of H 2 /N 2 - 1 A 1 combustion. Symbol key : + - i?2(*10), x - O 2 , o - H^O , 

A - H(* 2 • 10 5 ). 

3* Results 

The general behavior of the numerical simulations can be described as consisting 
of a brief initial period of chemical and fluid-dynamic adjustment, followed by an 
extended period of relaxation towards a perfectly mixed quiescent state. The simu- 
lations were run over a period of one to two initial turbulent timescales (r tj o), during 
which time the conserved scalar unmixedness was found to exponentially decay (see 
Fig. 3). During the same period the turbulent Reynolds number (determined using 
the mean viscosity at each time) decreased from ~ 60, to a value of ~ 40. 

Due to the chemical reactions taking place between the mixing fluids, the mean 
temperature and pressure typically rose by a factor of ~ 4/3 during the course 
of each simulation. From Fig. 3, the differences between the one- and two-step 
chemical calculations can be seen in terms of the mean species yield. It is evident 
that the one-step chemical description tends to underpredict the overall rate of 
reaction compared to the two-step case. Despite the fact that both chemical reaction 
mechanisms cause the system to tend to the same thermodynamic state in the 
absence of any mixing activity, it is apparent that the global reaction rate predicted 
by the one-step reduced mechanism is somewhat hindered in the presence of mixing. 

In both reaction mechanisms, the chain branching step H + O 2 ^ OH + O is 
the controlling component in the global rate. The one-step chemical mechanism 
determines the radical concentration from a quasi-steady state assumption that 
strictly only holds in the chemical equilibrium limit. Where mixing rates are high, 
fluid particles do not remain at the same equivalence ratio long enough to allow 
equilibrium conditions to be reached, and so the quasi-steady state assumption 
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FIGURE 3. Time history of unconditional mean statistics from uniform diffusivity 
simulations. (Left) Turbulent mixing quantities : + - ft, x - < e > (*2.5 ■ 10 10 ), o - 
| < u? > (*1 * 10 6 ). (Right) Major species mass fractions for each chemical case : 
+ - one step tf 2 , 0 - two step #2, n - one step 0 2 , 0 - two step 0 2 , x - one step 
H 2 0, A - two step H 2 0 

breaks down. The two-step mechanism on the other hand carries H as a computed 
species, and is not subject to this assumption. It would seem that the one-step 
prediction of H radical levels is lower than it should be under the mixing rates 
studied here, and this in turn limits the global reaction rate. 

The simulations with non-uniform molecular diffusivities displayed slightly differ- 
ent behavior to that discussed above. Discussion of non-uniform diffusivity effects 
is deferred to Section 3.2. 

3.1 Model results for uniform diffusivity cases 

The level of mixing intensity in the CMC model equations is described by the 
conditional mean scalar dissipation rate. It is important to accurately predict this 
quantity since it often closely balances the chemical production source terms (see 
Bilger 1989, Smith 1994). The conditional mean scalar dissipation rate profiles 
predicted by the model and observed in the companion simulations are plotted in 
Fig. 4 at three time stations t/rt = 1/3, 2/3, and 1. The conditional variance profiles 
of scalar dissipation rate are also plotted from the simulation data as an indication 
of the scatter in instantaneous dissipation rate from the conditional means. This 
scatter is not modeled and is a potential source of inaccuracy. 

It is apparent that despite the rather crude assumed-form PDF model used in 
computing the conditional mean scalar dissipation rate, the agreement with the 
observed profiles is quite reasonable. The most notable difference being the tendency 
of the predicted profiles to be greater than the observed profiles at very lean and 
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FIGURE 4. (Left) Comparison of modeled and observed conditional mean scalar 
dissipation rate <xli> profiles at different calculation times: + - DNS at t/r t = 
1/3, x - CMC at t/n = 1/3, o - DNS at t/r t = 2/3, A - CMC at t/r t = 2/3, n - DNS 
at t/rt = 1, and - CMC at t/T t ~ 1 (Right) Comparison of observed conditional 
mean and RMS profiles of scalar dissipation rate at two different calculation times: 
+ - mean at t/r t — 1/3, x - rms at t/T t = 1/3, o - mean at t/rt — 1, and A - rms 
at t/rt = 1. 

very rich mixture fractions. 

8.1.1 One-step chemical mechanism case 

In Fig. 5, conditional mean temperature and H 2 O mass fraction profiles are com- 
pared between CMC model predictions and simulation data. It is clear that the 
conditional mean profiles predicted by the model substantially exceed the profiles 
measured in the DNS. Further, the predicted profiles increase in magnitude with 
increasing time while the measured conditional mean profiles remain approximately 
stationary. 

The root mean square deviations from the conditional mean profiles, measured 
in the DNS, increase in magnitude with increasing time. These deviations axe not 
accounted for in the simple first order chemical closure currently used in the CMC 
model. It seems that these deviations from the conditional mean profiles are suffi- 
cient to cause the observed DNS conditional mean reaction rate to be substantially 
lower than the reaction rate modeled using conditional mean scalar profiles. 

Comparing the time histories of the one step chemistry model predictions and 
simulation results (see Fig. 6), it is apparent that the relative discrepancy between 
the unconditional mean profiles increases with time. The principal reason for this 
is that the mixing field slowly tends towards homogeneity, and thus conditional 
profile discrepancies near the mean mixture fraction become more prominent in the 
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FIGURE 5. Comparison of modeled and observed conditional statistics of temper- 
ature (left) and H 2 0 mass fraction (right) at different calculation times. + - DNS 
at t/rt — 1/3, x - CMC at t/rt = 1/3, o - DNS at t/r t = 1 , A - CMC at t/r t = 1 , 
n - DNS at t/r t = 1/3, and A - DNS at t/r t = 1 



FIGURE 6. Time history of unconditional mean chemical yields from one-step 
chemistry cases of model and simulation : (Left) + - DNS ^(*10), o - CMC 
H 2 {* 10), x - DNS H 2 0, A - CMC H 2 0, n - DNS 0 2 , 0 - CMC 0 2 . (Right) x - 
DNS pressure, A - CMC pressure, + - DNS temperature, o - CMC temperature. 



Modeling turbulent nonpremixed combustion 


311 



FIGURE 7. Comparison of modeled and observed conditional statistics for H radi- 
cal mass fraction at different calculation times. (Left) Mean modeled and observed 
profiles : + - DNS at t/r t = 1/3, x - CMC at t/r t = 1/3, o - DNS at t/r t = 1, A 
- CMC at t/rt = 1. (Right) Mean and RMS deviations from DNS : + - mean at 
t/rt = 1/3, x - mean at t/r t = 1, o - RMS dev. at t/r, = 1/3, A - RMS dev. at 

t/rt = 1. 


convolution with the PDF. 

Additionally there is a compounding effect of differences in mean pressure and 
temperature. The model tends to overpredict the heat release rate as a consequence 
of the first order chemical closure employed, and this in turn leads to overpredic- 
tions of temperature and mean pressure. This departure increases because the ele- 
vated temperatures and pressures cause even greater predicted heat release rates. It 
should be pointed out, however, that the model and simulation trends should con- 
verge given a sufficiently long time due to the limited amount of fuel and oxidizer 
present, and the fact that mixing motions will eventually disappear. 

3.1.2 Two-step chemical mechanism case 

The discrepancies between model and observation are significantly reduced in the 
cases considered with two-step chemistry. The major species profiles (not plotted) 
agree so closely as to be almost indistinguishable, save for the small perturbations 
associated with the DNS data. The only significant differences exist in the compar- 
ison of predictions and observations for the radical species (H ) and temperature. 

A comparison of conditional mean H mass fraction profiles for various calculation 
times is plotted on the left-hand side of Fig. 7. It is clear that the CMC model 
overpredicts the level of H present, but the relative degree of overprediction at the 
peak mass fraction decreases with time as the magnitude of the profiles decreases. 
It also seems that the overprediction of conditional mean scalar dissipation rate at 
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FIGURE 8. Comparison of modeled and observed conditional statistics for tem- 
perature at different calculation times. (Left) Modeled and observed mean profiles 
from two-step chemistry cases : + - DNS at t/r t = 1/3, o - DNS at f/r* = 1, x - 
CMC at t/r t = 1/3, A - CMC at t/r t = 1. (Right) Observed mean and RMS devi- 
ations at t/rt = 1 in one- and two-step chemistry cases : + - mean from two-step, 
x - mean from one-step, o RMS dev. from two-step, A - RMS dev. from one-step. 



FIGURE 9. Time history of unconditional mean chemical yields from model and 
simulation : n - DNS #(*10000), 0 - CMC #(*10000), o - DNS # 2 (*10), + - CMC 
# 2 (*10), x - DNS # 2 0, A - CMC # 2 0. 
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very rich mixture fractions, noted earlier, tends to cause the predicted H profile 
to be flattened at rich mixture fractions. In the right-hand plot of Fig. 7, the 
observed root mean square deviation profiles are plotted in comparison with the 
mean profiles. From this plot it appears that the absolute level of the deviational 
profiles decreases with time in accordance with the mean profiles. 

On the left-hand side of Fig. 8, conditional mean temperature profiles are com- 
pared for various times in the CMC and DNS two-step chemical calculations. The 
model profiles are somewhat higher than the DNS measured profiles; however, the 
difference is substantially less than that seen in the one-step chemistry case. Also 
in contrast to the one-step chemistry comparison, both sets of profiles increase in 
magnitude with time. This is another indication of the fact that the two-step for- 
mulation gives rise to a chemical system that is less perturbed by mixing processes 
when compared to a similar case with one-step chemistry. This fact is highlighted in 
the right-hand plot of Fig. 8, where conditional mean temperature data is compared 
between DNS simulations with one- and two-step chemistry at a time of t = r t . Not 
only is the mean profile greater in the two-step case, but the corresponding root 
mean square deviation profile is much lower. 

It appears that the more robust nature of the two-step chemical mechanism lends 
itself better to CMC modeling, in the cases studied here, than its one-step coun- 
terpart. This is because the two-step mechanism is less perturbed by the level of 
mixing intensity with consequently smaller mixing induced deviations from the con- 
ditionally averaged reactive scalar values. This reduction in the size of conditional 
deviations thus improves the accuracy of the conditional mean chemical closure. 

It is reasonable to assert that under more intense mixing conditions, that CMC 
models employing the two-step chemical mechanism would deviate to a larger degree 
from corresponding DNS observations. At very much higher mixing rates, the first 
order chemical closure would be invalidated altogether as the chemical system verges 
on extinction (see Bilger 1991, 1993). 

3.2 Observed differential diffusion behavior 

All of the DNS data and model predictions presented so far have been restricted 
to cases with uniform molecular diffusivity for all species (Le Q = 1, Pr = 0.75). 
The DNS data presented in this section was computed with constant non-uniform 
Lewis numbers determined from counterflow laminar diffusion flames (see Smooke 
1990). A Fickian diffusion approximation was used for all species, except nitrogen, 
which was the predominant background species for which the Lewis numbers were 
defined. 

One of the most notable aspects of comparing general simulation behavior, with 
and without differential diffusion, is the absence of a unique mixture fraction def- 
inition in the former case. Fig. 10 illustrates this fact by plotting the scatter of 
points for two different conserved scalars, one based upon the mass fraction of an 
inert species (iV 2 ) and the other based upon a combination of hydrogen and oxygen 
atomic mass fractions, at a time of t/rt = 1. Instead of adhering to the constant 
mixing line (unit-slope line passing through the origin), the computed points follow 
the characteristic sigmoidal (reversed in this case) trace of differentially diffused 
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FIGURE 10. Scatter plot of mixture fractions based on iV 2 mass fraction and H-0 
atomic mass fractions, for differential diffusion simulation at time t/r t = 1. 


conserved scalars. The reason for this behavior can be understood when it is re- 
vealed that the H — O mixture fraction has a positive linear relation to the light 
hydrogen bearing species and a negative linear relation to the heavier oxygen bear- 
ing species. Being lighter and more mobile than the other species, H and ff 2 tend 
to diffuse more rapidly to lean (lV 2 -based) mixture fractions than 0 2 can diffuse to 
rich (JV 2 - based) mixture fractions. The result is that the H — O mixture fraction 
values increases at lean iV 2 - based mixture fractions and simultaneously decrease at 
rich values. 

The chemical yields of the two step chemistry simulations are compared for cases 
with and without differential diffusion in Fig. 11. It is apparent that the differential 
diffusion case predicts a slightly greater reactant consumption rate, but with a 
less discernible increase in major production formation. It would appear that the 
additional reactants consumed by the differential diffusion case go towards creating 
the obvious excess of the radical species H . The temperature and pressure traces for 
the two simulations are very close; however, the differential diffusion case appears 
to have very slightly lower values. A clear difference is apparent between the two 
unmixedness traces (unmixedness of normalized iV 2 mass fraction), with the decay 
coefficient in the differential diffusion case being around ~ 0.86 of the coefficient in 
the uniform diffusivity case. 

Conditional statistics were calculated using an jV 2 -based mixture fraction defini- 
tion, and axe plotted in Figs. 12 and 13 for H radical mass fraction and temperature, 
respectively. It is apparent from Fig. 12 that the conditional mean H radical profiles 
from the differentially diffusive (diff-diff) case axe somewhat lower than the uniform 
diffusivity results on the rich side of stoichiometric. At the same time, the opposite 
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FIGURE 11. Time history of unconditional mean chemical yields from two-step 
chemistry simulations with and without differential diffusion. (Left) 0 - dd #(* 10 5 ), 
fl - ud H(* 10 5 ), o - dd # 2 (*10), + - ud #2(*10), A - dd # 2 0, x - ud # 2 0. (Right) 
0 - dd unmixedness, fl - ud unmixedness, o - dd temperature, + - ud temperature, 
A - dd pressure, x - ud pressure. 


is true on the lean side, where the differentially diffusive radicals have permeated 
this oxidizer-rich zone to a greater extent. 

The mean temperature profiles from the diff-diff cases reflect the greater incursion 
of H radical into the lean zone, in that they are significantly elevated over the 
uniform diffusivity profiles. The presence of greater radical numbers allows the 
exothermic global reaction to proceed at a more rapid rate, thereby liberating more 
heat. It is clear that the root mean square deviations in the diff-diff case are much 
greater than in the corresponding uniform diffusivity cases. Since temperature is 
strongly dependent on reaction activity and this is in turn dependent on radical 
availability, these temperature deviations may be related to the disparate mixing 
behavior of H radical and JV 2 (the conserved scalar) with the latter doing a poor 
job of tracking the mean transport of the former. 

Finally it is interesting to compare the iV 2 -based mixture fraction PDFs observed 
in each simulation case. The PDFs plotted in Fig. 14 are from the time t = r t and 
embody the main difference between the diff-diff and uniform diffusivity behavior. 
The diff-diff PDF equation has a non-zero source term arising out the definition of 
JV 2 mass fraction as the residual mass not accounted for by the reactive species mass 
fractions. This source term averages across mixture fraction space to provide a zero 
mean contribution, but serves to inflate the PDF at some mixture fractions. The 
source term (not plotted) has a sharp peak just to the lean side of stoichiometric, 
which results from H and # 2 incursion, and this peak causes the diff-diff case’s 
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FIGURE 12. Comparison of observed conditional mean (left) and RMS deviations 
(right) for H radical mass fraction in cases with and without differential diffusion 
at different calculation times : + * ud at t/r t = 1/3, x - dd at t/r t — 1/3, o - ud at 
t/r t — 1, A - dd at t/rt = 1. 



FIGURE 13. Comparison of observed conditional mean (left) and RMS deviations 
(right) for temperature in cases with and without differential diffusion at different 
calculation times : + - ud at t/r t = 1/3, x - dd at t/rt = 1/3, o - ud at t/rt = 1, 
A - dd at t/rt = 1. 
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FIGURE 14. Comparison of JV 2 - based mixture fraction PDFs from uniform and 
non-uniform molecular diffusivity simulations at time t = r t : + - uniform diffusivity, 
x - differential diffusivity 

PDF to increase at this mixture fraction. 

4. Discussion 

Much of current turbulent nonpremixed combustion modeling relies on the exis- 
tence of a unique mixture fraction, particularly the Joint PDF method (Pope 1985, 
1991, Chen et al. 1990), and, of course, the CMC method (Smith 1994, Smith et 
al. 1995). Rationales have been put forth that suggest differential diffusion effects 
are small at high Reynolds numbers; however, experimental evidence suggests that 
they axe substantial even in jet diffusion flames with Reynolds numbers as high as 
30000 (Smith et al. 1993). 

Let us briefly examine the impact of differential diffusion on the CMC method and 
discuss the additional modeling issues which arise. The equation for a composite 
conserved scalar (£ = X)a=i a «^o) is given by, 

+ = + (15) 

where is a differential diffusion source term equal to 

MO. -0 { )§|). (16) 

J a— l J 

From this equation it can be shown that the corresponding PDF equation is given 

by, 


d d 

t^(< p\rj> P„) + — (< p Ui | T) > P v ) + et 


( 17 ) 
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and following the methodology of Klimenko (1990) yields the following CMC equa- 
tion. 


<p\rj> ^f-+ < P u i I V > = \<PX\V> ^jr+ <pw a \r,> +e* + e* Q 


The residual terms e* and e*q axe defined as, 


(IS) 


e; = i v > p,) - §j(< *«*'. i -? > p,) i <i») 

«r <20) 

The additional terms in the CMC and PDF equations are dependent on the defi- 
nition of mixture fraction selected els the conditioning variable. Note that the eg 
term becomes small at high Reynolds numbers, but it is not clear that the same 
can be said for e*. 

The choice of a conserved scalar els a conditioning variable is governed by two 
criteria. Firstly, the scalar should be representative of the molecular transport of as 
many importsmt chemical species els possible so as to minimize deviations from the 
resultant conditional averages. Secondly, in order to be able to model the evolution 
of the conserved scalar PDF (and thereby determine < px I V >) m a simple manner, 
the mixture fraction source term if should be as small els possible. These criteria 
may prove to be conflicting. For example, it is possible to use an inert tracer species 
as a conserved scalar, thereby malcing if identically zero, but as this definition does 
not include any of the reactive species that are being tracked, deviations from the 
conditional means may be too great to effect a chemical source term closure. 

Klimenko (1994) provides an equation for the conditional mean square deviation 
6 Q from a conditional mean reactive scalar Q a in isotropic turbulence. This equation 
is slightly modified in the presence of differential molecular diffusivity to become, 


. ee Q i . de Q ., , 
<p ^ > ~dT~ 2 <PXIV> 'W =< pWaVa 


t] > — 2 < D a {-^-) 2 | r] > +e£+e@. 

( 21 ) 


where the residual terms are analogous to e* and eg but involve the conditional 
mean deviation rather than the conditional mean. The instantaneous change and 
transport of the conditional mean deviation is balanced against a chemical-instability 
source term, a deviational dissipation term, and residuals. Where the Reynolds 
number is moderate and the choice of conserved scalar is poor, the residuals will 
tend to increase the level of deviations. When combined with the nonlinear ampli- 
fication provided by the chemical term, this added source of conditional deviation 
can cause levels to increase substEmtially (see Section 3.2) Emd thereby invalidate 
any first order chemical closure. 

It may be that the chemical source term < pw Q \ 77 > can be closed using a second 
order method such as that applied by Li and Bilger (1993) to atmospheric pollutant 
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reactions in a turbulent mixing layer. However, in that case the chemistry was 
isothermal and one-step in nature and did not have significant differential diffusion. 
At higher Reynolds numbers the effect of differential diffusion is diminished, and it 
may be that practical devices that can be highly turbulent do not require a diff- 
diff treatment. At higher Reynolds numbers, however, deviations arise in Eq. 21 
because of mixing interference via the chemical instability term (see Section 3.1), 
and it may be necessary to develop doubly conditional moment closures in order to 
model these conditions. 

5* Conclusions and future plans 

In this study, predictions from the CMC method for modeling turbulent non- 
premixed combustion were compared to DNS data for hydrogen burning in an 
isotropic decaying turbulent field. One- and two-step chemical mechanisms were 
used in both the model and simulation in order to study the effect of chemical 
complexity upon first order CMC chemical closure. 

It was found that the one-step chemical mechanism was hindered to a greater 
extent over the two-step mechanism, under identical mixing conditions, as a result 
of a breakdown in the one-step assumption for radical partial equilibrium. This 
interference by the mixing processes lead to larger deviations from conditional mean 
reactive scalar profiles. This in turn made the one-step chemical system harder to 
model with the CMC method than the two-step system under the same mixing 
conditions. 

The addition of differential molecular diffusivity to the analysis tended to increase 
the level of reactive scalar deviations from conditional means under the conditions 
studied. The different rates of species transport tended to modify the overall rate of 
chemical reaction in the hydrogen system. The lack of a unique conserved scalar as 
a conditioning variable caused conditional mean scalar profiles to shift in mixture 
fraction space according to the choice of conserved scalar. 

It was suggested that the increase in conditional deviations that arise from differ- 
ential diffusion effects is a potential source of serious in implementing a conditionally 
averaged first order chemical closure. 

The future plan for this modeling project can be outlined as follows: 

• Perform three-dimensional simulation under same kind of conditions to include 
the vortex stretching mechanism and improve the size of the statistical sample. 

• Develop model refinements for treating differential diffusion — may require solv- 
ing for conditional deviations. 

• Investigate doubly- conditional moment closure methods in spatially degenerate 
case. This will allow high intensity near-extinction behavior to be examined. 

• Use flexible chemical module for other mechanisms such as H<i — CO. 
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